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SECTION  1 


INTRODUCTION 


The  hyperbolic  nature  of  the  Euler  equations  in  supersonic  flow  permits  a 
solution  by  spatial  marching.  Here  an  initially  known  flow  field  on  a 
crossflow  plane  near  the  missile  nose  is  marched  down  the  missile  length.  For 
a  sharp  tipped  missile  the  initial  flow  field  can  be  constructed  from  a 
conical  solution  while  a  blunt  body  solver  must  be  used  in  the  case  of  a  blunt 
shape.  The  computation  terminates  when  the  base  of  the  missile  is  reached  or 
when  subsonic  regions  are  encountered. 

This  report  describes  the  application  of  a  second  order  Godunov  method  to 
tactical  missile  shapes  (i.e. ,  bodies  with  thin,  low  aspect  ratio  lifting 
surfaces).  The  grid  is  generated  using  a  multiple  zone  method  that  divides 
the  crossflow  plane  up  into  several  quadrilateral  zones.  The  body,  fin,  and 
tracked  bow  shock  surfaces  must  lie  on  the  zone  edges.  A  typical  zone 
structure  for  a  body-wing  missile  is  illustrated  in  Figure  1  while  the 
quadri lateral  zone  geometry  is  shown  in  Figure  2.  The  numerical  method  is 
cast  in  control  volume  form  and  consists  of  a  predictor  and  corrector  step. 

The  predictor  step  advances  the  primitive  variables  using  Euler's  equations  in 
non-conservation  form.  Derivatives  are  computed  using  a  limited  central 
differencing  procedure.  The  corrector  step  modifies  Godunov's  method  by 
assuming  linear  property  variations  within  each  control  volume. 

The  multiple  zone  computational  procedure  has  been  previously  applied  to 
missiles  using  a  MacCormack  type  scheme  with  characteri sti c  boundary 
conditions.1’^  At  wing  edges  and  surface  discontinuities,  these  boundary 
conditions  are  appended  with  special  procedures  which  apply  oblique  shocks  or 
expansions  as  appropriate.  These  methods  have  been  successfully  applied  to 
missiles,  but  lack  robustness.  In  simple  cases  (e.g.,  a  tangent  ogive  at  low 
incidence),  accurate,  efficient  solutions  are  achieved.  However,  for  more 
complicated  shapes  the  user  must  adjust  the  artificial  viscosity  levels  and 
special  procedures.  A  first  order  Godunov  method  has  also  been  developed 
using  the  multiple  zone  approach’  which  is  extremely  robust  and  devoid  of 
special  procedures.  In  regions  where  shocks  occur  these  results  are 
comparable  to  those  obtained  by  the  MacCormack  scheme.  However,  on  smooth 
bodies,  accurate  prediction  of  surface  pressures  can  only  be  obtained  with 
extremely  fine  meshes.  This  general  recipe  for  constructing  a  second  order 
accurate  Godunov  scheme  was  suggested  in  Reference  4.  It  has  been  applied  to 
unsteady,  compressible  flow  in  References  5,  6  and  7  and  to  steady  supersonic 
flow  in  References  3  and  9. 

The  second  order  Godunov  scheme  presented  in  this  paper  is  Imbedded  in  a 
multiple  zone  framework  which  features  the  division  of  the  computational 
domain  into  quadrilateral  zones  of  the  type  shown  in  Figure  2.  The  edges  of 


adjacent  zones  spacially  coincide  with  one  another  and  do  not  overlap.  The 
quadri 1 ateral  zone  must  be  consistent  with  the  following  constraints  (see 
Figure  3  for  zone  numbering  convention): 

1.  Edge  2  can  only  abut  to  edge  4.  Edges  1  and  3  cannot  abut  to 
other  edges. 

2.  Surfaces  and  symmetry  planes  must  lie  on  zone  edges. 

3.  Bow  shocks  which  are  to  be  fitted  must  lie  on  edge  3. 

The  geometry  of  the  zone  edges  can  be  defined  using  either  cylindrical  or 
cartesian  coordinates.  These  constraints  are  sufficiently  loose  to  allow  a 
large  number  of  different  internal  or  external  flows  to  be  handled.  This 
includes  body-wing-tail  missiles,  wings  alone  and. inlet  flows.  Geometries 
including  features  such  as  detached  inlets  or  a  tail  located  on  a  wing  require 
abutting  edges  1  and  3,  as  well  as  2  and  4.  These  types  of  shapes  cannot  be 
handled  by  the  present  second  order  Godunov  code.  Reference  2  describes  a 
MacCormack  based  multiple  zone  method  for  computing  such  flows. 

A  user  guide  for  the  method  described  in  this  report  is  provided  in 
Reference  9. 
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SECTION  2 

THE  RIEMANN  PROBLEM 


In  steady  supersonic  flow  the  Riemann  problem  represents  the  confluence 
of  two,  two-dimensional  supersonic  streams,  as  is  illustrated  in  Figure  4.  At 
tne  point  of  stream  intersection,  shocks  and/or  expansion  fans  form,  turning 
both  streams  to  a  common  direction.  The  appropriate  direction  is  the  one 
producing  the  same  pressure  in  both  streams.  The  two  final  streams  need  not 
feature  the  same  density  or  velocity  and  a  slip  line  usually  forms  between 
them.  Solutions  are  self-similar  in  z  and  feature  constant  properties  along 
any  line  passing  through  the  point  of  initial  stream  intersection. 

Figure  3  depicts  the  Riemann  problem  formed  by  the  intersection  of  two 
streams  with  properties  (p+,  p  +,  u+,  w+)  and  (p-,  p-,  u-,  w-).  Upward 
deflection  of  the  (+)  stream  is  described  by  the  shock  relations  and  leads  to 
an  increase  in  pressure  while  downward  deflection  is  associated  with  an 
expansion  and  leads  to  a  decrease  in  pressure.  On  the  other  hand,  upward 
deflection  of  the  (-)  stream  leads  to  an  expansion  and  downward  deflection  to 
a  shock.  To  solve  the  Riemann  problem,  it  is  necessary  to  determine  the  final 
flow  direction,  which  produces  the  same  final  pressure  in  both  streams.  The 
shock  and  expansion  equations  define  the  relation  between  the  two  initial  flow 
directions,  6  and  pressure,  p  ,  and  final  pressure  and  direction,  p,  5. 

These  can  be  written  expressing  5  as  a  function  of  <$  ,  p  ,  and  p: 


5f  =  5+  ±  tan' 


(Pf-p±)  l2  rp±(27M±2-nl)  -  (Y+I)pf-|1/2 


(Pf-P±j  -j  |- 

y^p±"pfj  . 


=  <5+  t  v(M+)  *  v(Mf ) 


1y+I)pf+TY-l 


:expansi on 
pf  <  P+ 


■]  )■ 


shock  (1) 
Pf  >  P+ 


T-l 
.p*.  y 


.here:  Mf  .  (1  f  )  -  1] 

u(M)  -  t  -  ta„-V-l]1/2 

C2 

c  - 


* 


*V*7 


mm'w 


Here  the  upper  and  lower  signs  are  associated  with  the  (+)  and  (-)  stream 
respectively.  The  resulting  (p  -  5)  curves  for  the  +  and  -  streams  are 
graphed  in  Figure  4.  The  solution  to  the  Riemann  problem  corresponds  to  the 
intersection  of  the  two  illustrated  curves. 


2.1  COMPLETE  SOLUTION 

The  non-1 i nearity  of  Eqs.  (1)  -  (2)  precludes  a  closed  form  solution. 
Instead,  an  iterative  procedure  is  used  which  proceeds  from  step  n  to  n+1  as 
f  ol  1  ows: 


1.  At  the  start  of  iteration  step  n+1,  two  points  are  known  on  both 

the  +  and  -  (p-<5)  curves:  (pj,  fi”),  (p""1,  s"'1).  The  equations  for 

the  lines  passing  through  the  two  known  points  on  each  curve  are 
determi ned. 

2.  The  intersection  of  the  two  lines  is  determined  and  p  is  taken  to 
be  the  pressure  at  the  point  of  intersection. 

3.  and  <5_  are  determined  by  evaluating  Eqs.  (1)  and  (2)  using  the 
value  of  p  determined  in  step  2. 

4.  If  j  <5+-<5_J  <  10~3,  the  iteration  is  terminated. 


To  start  the  iteration: 


(3) 


Here  p|_  is  the  linearized  pressure  solution  discussed  in  the  next  section. 

<5^  are  determined  by  evaluating  Eqs.  (1)  and  (2)  using  P|_. 

The  Godunov  procedure,  which  is  described  in  Section  3.2,  requires 
knowledge  of  the  flow  properties  associated  with  a  specific  direction,  9  .  As 
indicated  in  Figure  3,  three  different  regions,  (R±,  F,  E)  are  featured  on  the 
expansion 

side  of  a  slip  line  and  two  (R+,  S)  on  the  shock  side.  If  9  lies  in  R+,  the 

initial  conditions  provide  the  needed  properties.  For  9  lying  in  S  or  E,  the 
variables  p,  5  are  determined  by  the  Riemann  solution  and  the  density  is 
computed  using  the  oblique  shock  relations: 


( y+1 )M2si  n29$ 

™— A  "M  ' 


P 


(4) 
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fi&S 

O  * 
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where:  sin29s  =  [(x+1)  ?  +  (t-1)]/ (2fM2) 


(5) 


In  region  E,  the  isentropic  relations  provide  the  correct  density.  However, 
the  shock  relation  is  less  expensive  to  apply  and  yields  the  required 
accuracy.  When  8  lies  in  F,  the  flow  direction  6  along  0  is  given  by: 


5  =  0  ?  si n-1  {-L.} 


(6) 


Here  M*  is  the  Mach  number  along  the  direction  8.  Using.  Eq.  (2): 


5.  t  v(M  )  =  9  ±  v(M* )  +  sin'1^} 
*  1  M 


and  solving  for  M  yields: 


(7) 

(8) 


*  —o  7  +  v(M  )  +  (5-0)  y/2 

M  =[1+Ftan2{ - ±- - ± - }]1/2 


With  M  known,  the  remaining  properties  along  9  follow  immediately  from  the 
isentropic  relations. 


The  angle  between  the  shock  and  the  slip  line,  9  ,  is  given  in  Eq.  (5), 

while  the  angle  between  the  slip  line  and  the  head  ana  tail  characteri sti cs  of 
the  expansion  fan  are  the  characteristic  angles,  based  on  the  Mach  number 
upstream  and  downstream  of  the  fan  respectively. 


2.2  APPROXIMATE  SOLUTIONS 


Approximate,  closed  form  solutions  to  the  Riemann  problem  can  be  obtained 
in  cases  featuring  similar  +  and  -  states.  This  often  circumvents  the  need  to 
use  an  iterative  method.  To  develop  approximate  solutions,  the  shock 
relation,  Eq.  (1),  is  expanded  in  the  following  Taylor  series: 


P  =  P±±k.  (5-5  )  +  k  (5-5+)2  i  k  (5-5+)3  +  0[(5-5j4] 


(9) 


where:  ^  =  ptYM2/ (M2-!) 1/2 


k2  = 


(nl)M4  -  4(M2- 1 ) 


4(M^-1  )~2 
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A 


V'. 


(M^-l) 7/2 


(y+1)'  m8 

^rr~  ± 


(7  +  12y  -3/)  w6 

- n - ± 


+  5  (y*1)mJ  -  mJ  +  | 


The  Prandtl -Meyer  expansion  series  is  identical  to  the  above  through  second 
order.  The  third  order  terms,  for  y=1.4,  are  typically  within  10%  of  each 
other,  and  Eq.  (9)  will  be  used  to  represent  both  shocks  and  expansions. 
Retaining  only  linear  terms  and  requiring  that  the  final  upper  and  lower 
stream  pressures  be  equal  yields  the  linear  solution: 


PL  -  P+  + 


><l+{[P_-p+]  +  [<5_-<5+]} 


I<!  +  <i  ) 


1 


(10) 


((P.“P+)  +  ICj  (5_-6+)} 


6L  =  5+  + 


+  K.J 

+ 


This  estimate  can  be  improved  by  considering  the  second  order  term  in  Eq. 


;  oj 


La* 


be  the  quadratic  solution  and  write: 


5g  “  5+  =  (5q  ~  *5^)  +  ($[_  -  <5+)  =  A6q  +  A6l 


(ID 


Substituting  the  above  into  Eq.  (9),  and  dropping  the  ASq  term  results  in: 


p  =  PL  +  k2  a^_£  +  A5QC±k1  +  2k2  A6L  ] 


(12) 


Equating  the  upper  and  lower  stream  pressures  yields: 


2 


A6Q  =  (k2  a^  c  -  k2  a^  )/ (2k2  a^  -  2k2  a^  +  kj  +  kj  ) 


(13) 


•+  ~+ 

The  final  pressure  can  be  evaluated  to  second  or  third  order  using  Eq.  9. 


Properties  along  a  direction  9  are  evaluated  using  the  same  procedure 
discussed  in  the  preceding  section,  except  in  the  linear  case.  Here  the 
Riemann  solution  is  assumed  to  contain  two  distinct  regions  separated  by  lines 
along  which  flow  properties  are  discontinuous,  as  is  illustrated  in  Figure 
5.  In  particular,  it  is  assumed  that  expansion  fans  are  infinitely  thin. 


The  efficiency  of  the  Riemann  problem  strongly  impacts  that  of  the 
Godunov  procedure.  To  minimize  the  computational  resources  devoted  to  the 
Riemann  problem,  the  approximate  solution  is  calculated  as  follows: 


J Cru 


ftAkyrj-'Vj 
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1.  The  linear  pressure  and  flow  direction  are  determined  from  Eq.  7. 

2.  If  max  (j<5^-<5+),  j<SL-<sJ)  <  *02,  the  linear  solution  is  used  as  the 
final  Riemann  solution. 


3. 


4. 


5. 


If  max  ( |  <S|_-<5+| ,  |  | )  >  .02,  the  quadratic  solution  is  evaluated 

from  Eq.  13. 


If  .02  <  max  (|5q-6+|,  |  5q-sJ  )  <  .06,  the  pressure,  pq,  is  evaluated 

using  first  and  second  order  terms  of  Eq.  9.  Slightly  different  pq 
values  are  obtained  from  the  upper  and  lower  streams.  The  final 
pressure  is  an  average  of  pq  . 


If  .06  <  max  (|$q-$+|,  |5q-sJ)  <  .10,  the  first,  second  and  third 
order  terms  of  Eq.  9  are  used  to  calculate  pressure,  dc.  pc+  and  pc 


are  averaged  to  determine  the  final  pressure. 


6.  If  max  (|5q-5+|,  | 6  - 6_ | )  >  .10,  the  iterative  procedure  outlined  in 
the  last  section  is  applied. 

2.3  APPROXIMATE  PROBLEM 

An  alternative  to  obtaining  an  approximate  solution  to  the  supersonic 
Riemann  problem  is  to  define  an  approximate  Riemann  problem  which  can  be 
readily  solved.  Following  the  general  approach  outlined  by  Davis  for 
conservation  laws ,  a  simplified  Riemann  problem  for  supersonic  flow  can  be 
constructed.'  It  features  a  single  set  of  intermediate  properties  separating 

the  two  initial  states  as  is  shown  in  Figure  6.  The  slopes,  a+  and  a',  of  the 
two  surfaces  across  which  flow  properties  are  discontinuous,  are  defined  as 
follows: 

a+  =  min(A_+,  A_“) ,  a”  *  max(\++,  \+~)  (14) 


where:  the  subscript  refers- to  the  root 

the  superscript  refers  to  the  +  or  -  stream 


wu+a 


aVu2+w2 


7 


(15) 


The  parameters  a±  are  easily  computed  and  ensure  that  the  complete  Riemann 

problem  wave  speeds  are  faster  than  a_  and  slower  than  a+.  This  assures  that 
the  entropy  condition  is  satisfied  and  that  the  method  converges  to  the 
correct  physical  solution. 


The  intermediate  state,  denoted  by  a  superscript  asterisk,  is  calculated 
by  balancing  mass  and  momentum  fluxes  through  the  dashed  control  volumes  of 
Figure  6.  A  quadratic  equation  must  be  solved  to  determine  the  primitive 
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variables  p,  o,  u,  w.  However,  the  information  necessary  for  the  Godunov 
scheme  is  simply  the  fluxes  normal  to  some  direction,  9,  as  shown  in  Figure  6: 


★  ★ 
p  vn 


*  *  +  # 

P  V  w  +p  cos  0 
n 


(16) 


L*  *  *  * 

p  v  u  +p  sine 

n  J 

Here  vn*  is  the  velocity  normal  to  the  9  =  constant  plane. 


— * 

p  can  be  computed 


directly  without  the  need  to  determine  thd"  primitive  variables.  To  accomplish 
this,  conservation  is  satisfied  for  the  two  dotted  control  volumes  shown  in 
Figure  7.  This  yields  the  following  vector  equation: 


a+U  =  aV  - 


cos  9  cos"§ 


.  *  -  -  F  F 

-a  u  «  -a  U  +  — s 

COS  0  COS  9 


where:  a+  *  a+  -  tane 

a“  »  a"  -  tane 
Solving  for  F  produces: 


upper  control  volume 
lower  control  volume 


(17) 


F*  = 


a'a'^U'1'  -  U,  )cos9  +  a^F* 
(a+  -  a") 


(18) 


When  applying  the  approximate  Riemann  solution  to  the  Godunov  scheme,  it 
is  necessary  to  determine  F  along  some  arbitrary  direction  9. 

If  a"  <  tan9  <  a+,  the  above  flux  definition  is  the  appropriate  one.  However, 

for  tane  >  a+  or  tane  <  a”,  the  F+  and  F“  fluxes,  respectively,  must  be 
used.  This  can  be  accomplished  using  the  above  formula  if  a+  and  a"  are 
redefined  as  follows: 


a+  *  max(a+  -  tane,  0) 
a"  *  min(a’  -  tane,  0) 


(19) 
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SECTION  3 


NUMERICAL  METHOD 


3.1  ZONE  DESCRIPTION 

In  the  crossflow  planes  z  =  constant,  zones  are  generalized 
quadrilaterals  as  shown  in  Figure  2.  Each  zone  is  described  with  respect 
to  (s,  t,  v)  coordinates,  which  either  represent  cylindrical  (r,  <t>,  z)  or 
cartesian  (x,  y,  z)  coordinates  (see  Figure  7).  Cylindrical  coordinates 
facilitate  treatment  of  wing-body  configurations,  while  cartesian  coordinates 
can  be  applied  to  wing  alone  cases.  The  numbering  system  used  to  designate 
the  edges  and  corners  of  a  zone  are  indicated  in  Figure  2.  The  locations  of 
edges  1  and  3  are  defined  by  the  functions  b(x,  v),  while  the  coordinates  of 
edges  2  and  4  are  defined  by  '?( s ,  v),  and  a(s,  v),  respectively. 

When  (s,  i,  v)  represent  cylindrical  coordinates,  the  bow  shock  may  be  fitted, 
but  it  must  be  located  on  edge  3. 

The  mesh  in  each  zone  is  defined  by  the  transformation  T 1  •  T2,  where  T1-: 
C5.  n,  0  -  (s,  t,  v)  and  T2  (s,  t,  v)  (x,  y,  z).  T2  is  given  by: 

x  =  s  cos  (t) 
y  =  s  sin  (t) 
z  =  v 

x  =  s 
y  =  x 
z  =  v 

while  T]_  is  defined  as: 

s  *  b(  t*  ,  c) 
x  =  a( s' ,  ?) 
v  »  4 

where 

t1  *  t4U)  +  (^(c)  -  x4U))g(n) 

T"  3  t3(c)  +  (t2(j ;)  -  T3(c))g(n) 

s’  =  s4(c)  +  (s3U)  -  s4(c))f(5) 


cylindrical  coordinates 


cartesian  coordinates 


(20) 


+  Cc(t",  c) 

+  C*(s",  0 


b( t‘  ,  ?)]f(c) 
o( s',  c)Jg(n) 


(21) 


and  (s.j ,  -r.j),  i  =  1,  2,  3,  4  are  the  coordinates  of  the  corners.  Here  f  and  g 

are  mesh  clustering  functions  with  f(Q)  -•  g(Q)  =  0  and  f(l)  =  g(l)  =  1.  The 
computational  domain  for  each  zone  is  bounded  by  0  <  g  <  1,  0  <  n  <  1,  and 
?  >  0. 

In  each  crossflow  plane,  the  metric  quantities,  g  ,  g  ,  g  ,  ny,nu,  n_, 

x  y  4  x  y  4 

must  be  evaluated  at  cell  centers.  These  parameters  are  used  in  the  predictor 
step  and  to  compute  the  step  size.  This  is  accomplished  analytically  as 

fol 1 ows: 


(Vx  *  Vx)/j  ;  nx 
(Tnsy  "  Vy)/j  ;  ny 
(snT?  '  V?)/j  ;  nz 


("Tgsx  +  sgTx)/j  ;  ?x  =  ux 
("Vy  +  sgT y)j  ;  cy  =  \ 
(V?'  Vs)/j  ;  cz  3  vz 


where:  j  =  s5Tn  “  snT5 

The  x,  y,  z  derivatives  of  s,  t,  v  are  given  in  cylindrical  and  cartesian 
coordinates  by: 


cos(t),  sy  *  sin(t),  sz  *  0 
cos (  t) /s  ,  t  =  -si  n(  t)/s  ,  T 
^•0.vx.l 


cylindrical  coordinates  (23a) 


i,  sy  =  sz  =  0 

1  ,  Tx  =  Tz  •  0 

c  ux  *  uy  *  0 


cartesian  coordinates 


(23b) 


The  derivatives  of  s,  t,  with  respect  to  g,  n,  are  ■eval uated 
analytically.  The  expressions  for  these  quantities  involve  corner  values  of 
sp  and  xp,  which  can  be  computed  using  the  formulas  given  below: 

S1  =  (*-bT*s);  =  ( V*Sbv)/(l“bT*s)  corner  1 


(c  v+ct4»v)/  ( 1  -c  t4»s  ) ;  t2  =  ( ^v+^scv)/ (l-cTtps) 


corner  2 


*  <cv  +  ctaJ/ ^"c-ras^  ;  t3  3  K  +  ascv^^1_CTCTs^  corner  3 
3  (bv  +  bTav)/(1  -bTffs)  ;  t4  ■  (bv  +  asb v)/ ( I  -  bTas)  corner  4 
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Using  the  corner  values  of  s  and  t  ,  derivatives  of  s',  s",  x‘  and  t"  can  be 
computed  as  follows:  ?  ? 


s' 5  ■  ( s 3  -  s4)f ' (C) 

s' ,  =  s4  7(C)  +  s3  f ( 0 

4  c  c 

s'  =0 
n 

t'  =0 
C 

t' ?  =  t4  g(n)  +  t1  g(n) 

t'„  -  <ti  -  Mi?'  <'> 


s”5  ■  (s2  -  Sl)f  (5) 

s"c  =  s1  7(C)  +  s2  f(c) 
c  c 

s"  =  0 
n 

t"  k  =  0 

t‘‘  =  g(n)  +  t?  g(n) 
T” n  3  ( t2  ’  T3)g'(n) 


where: 


f'U)  ,  g'(5)  -  Mil  ,  f(S)  -  1-fU),  gU)  *  l-g(C) 


1  —  »  y  1  — ^ —  »  '  \ s /  =  x-1  »  yv^  -  i-y^/ 

The  general  formulas  for  the  transformation  quantities  then  follow: 

s^  =  (c-b)f' (C) 

sn  =  vV(^)  +  vY<«> 

s  =  (b  t'  +  b  )f(c)  +  (c  t"  ♦  c  )f(c)  ' 
c  t  z  v  t  ?  v  (26) 

Tc  =  ars'  ^(n)  +  ^s^ ?  3^ 
s  (*-o)g‘ (n) 

t5  -  (v‘t  +  +  ^ss"t  +  ♦vW*) 

Additional  mesh-related  information  required  by  the  Godunov  procedure  is 
the  area  of  each  cell  in  the  crossflow  plane  (A^),  and  vectors  normal  (n)  to 
the  sides  of  the  control  volume  (see  Figure  8).  JThese  quantities  are  computed 
as  follows: 


Aij  *  7  V64  x  V51 


'i  +1/2 ,j  '  7  v*24 


(V« -  x  V,« ) 


(27a) 

(27b) 


Here,  V,,  =  (x.  -  x  ,  y.  -  y,,  z^  -  zO,  where  the  subscripts  refer  to  the 

J  *  J  1  *  J 

corner  numbers,  which  for  cell  (i,j)  and  edge  (i  +  1/2,  j),  are  defined  in 
Figure  8.  If  the  edge  corners  are  not  in  a  plane,  n  lies  in  an  area  weighted 
average  normal  direction  to  the  edge  (see  Reference  12).  Note  that  the 
magnitude  of  the  normal  vector  is  the  side  area. 


'A  V 


•  wy 

■  ,  -1  .  .  *>.  ,*'*  -•*  S  v  . *•.-  V  -i  "i  ■ 


vWV'' ‘V* v.f 
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3.2  SECOND  ORDER  GODUNOV  SCHEME 


Using  the  notation  and  coordinates  of  Figure  8,  control  volume  mass  and 
momentum  conservation  equations  are  given  by: 


TfTl+1  |  .11  p  p  p  T 

i,j  =  Ui,j  ‘  h  i  -*-1/2  ,j  +  hi  -1/2, j  "  hi  ,j+l/2  +  Fi,j-l/2 


0n  .  =  An  . 


g 

r  -i  i 

"  pV 

PW  +  P 

u2 

pwV  +  n2p 

pwu 

= 

u3 

Fi  +1/2  ,j  “ 

puV  +  nYp 

PWV 

u4 

pvV  +  up 

- 

i  J 

i ,  j 

«  * 

V  =  ni+l/2,j(u’  V*  w)i+l/2.j 


Here  U  is  the  flux  in  the  z  direction  which  passes  through  the  shaded  cell 
ends  (see  Figure  8)  while  the  F’s  are  the  fluxes  associated  with  the  remaining 
cell  edges.  Eqs.  (28)  are  closed  using  the  constant  total  enthalpy  condition 
and  the  perfect  gas  equation  of  state  which  yield  the  constraint: 


(28) 


->  i  +  l/2,j 


BpT ,Vi<«2*v2 

In  non-conservation  form,  Eqs.  (28)  becomes: 

u  (Sx»Vnxpn?  1 

?  VV  [  P  J 


(29) 


_-lFvSU  +  \V  +  VlVnj 


(30) 


■g-^-g-  LH^  -  wH^  +  P/p] 


p  =*  — v-r-  C-pwH,  +  pa2H2  -  wP]  -  5  p  -  u  p 
?  (w2-^)  1  2  z  5  z  n 


U(V?  +  Vn}  +  V(SP?  +  Vn) 


I 


y 

k 
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H.  *  a  {£  u_  +  £  v_  +  5  w,  +  n  u  +  n  v  +  n_w  } 
l  x  £  y  ?  z<i  xn  yn  z  n 

H0  =  Dw_  +  Vw 

The  solution  to  a  supersonic  flow  field  can  be  marched  in  the  z  direction 
using  Eq.  (28).  The  complete  second  order  Godunov  method  evaluates  the  F 
fluxes  and  completes  a  marching  step  as  follows: 

1.  Derivatives  of  p,  p,  u,  v,  w  are  computed  using  a  limiter.  To 

illustrate  the  limiter,  consider  the  derivative  which  is  calculated 
as  follows: 


if  (pi  +  l,j  "  pi,j^pi,j  ”  pi  -l,j  ^  < 


3p, 

TV  i ,  j 


(F)tni 

A£ 


n  lpi+l,j 


otherwi se 


±ML,  K  . j  I  ’  Klpi,j 


where:  F  «  sign(p1+1_j  -  pf  j) 


Here  K  is  an  adjustable  constant,  which  is  normally  set  to  unity  at 
interior  points. 

The  metrics  z,  ,  £  ,  £  ,  are  calculated  at  each  cell  center  in  the 
x  y  z 

computation  using  Eqs.  22. 

The  step  size,  Az,  is  determined  from  the  (CFL)  condition  which  is 
described  below  and  discussed  in  Appendix  A. 

The  predictor  values,  p,  p,  u,  v,  w  are  calculated  at  z  +  az/ 2  using 
the  derivatives  determined  in  step  1  and  the  metrics  from  step  2 
applied  to  Eq.  (30). 

The  coordinates  of  the  control  volume  corners  are  determined  from 
Eqs.  (21),  and  the  vectors  normal  to  cell  edges  are  computed  using 
Eq.  (27b). 

Properties  adjacent  to  cell  edges  at  z  +  Az/2  are  calculated  using 
the  cell  center  predicted  values  from  step  4  and  the  derivatives 
from  step  1.  This  produces  two  sets  of  properties  associated  with 
each  cell  edge;  one  for  each  of  the  adjacent  cells.  For  example, 
the  two  pressures  on  the  lower  edges  of  cell  (i,  j)  in  Figure  8  are: 


pt,j  -  7 -ft  i,j  from  cel  1  (i  ,  j ) 


H-t  .i  +  \  lei  i-l.i 


from  cell  ( i - 1 ,  j ) 


•ftj 

‘Vn 
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A  Riemann  problem  is  constructed  at  each  cell  edge  using  the  two 
sets  of  properties  determined  in  step  6.  The  two  dimensions 
associated  with  the  Riemann  problem  are  defined  by  plane  R  which 
contains  the  cell  edge  normal  vector  and  is  aligned  with  the  2 
axis.  The  two  initial  states  consist  of  (p,  ,  vn,  w) ,  at  the  two 

adjacent  cells,  where  vp  =  (nx,  n  )  •  (u,  v)/j(nx,  ny)|.  Also 

computed  for  each  initial  state  is 

vt  =  (“nv»  nx)  *  (u»  v)/|(n  ,  nil  which  is  the  velocity  component 
tangent  to  the  cell  edge  trace  in  the  z  =  z  plane. 

The  Riemann  problem  associated  with  each  cell  edge  is  solved  using 
the  techniques  described  in  Section  2.2. 

The  angular  orientation  of  the  cell  edge  in  plane  R  is  calculated 
with  respect  to  edge  (i+l/2,i).  This  orientation  is: 


Cell  edge  properties  are  constructed  using  the  Riemann  solution 
along  direction  9,  and  by  specifying  vt.  In  terms  of  the  properties 
provided  by  the  Riemann  solution,  p0,  o0,  q0,  60,  the  cell  edge 

properties,  denoted  by  a  subscript  e,  are  given  by: 


ue  =  (q0sin<5nx  -  ny ) / 1 nx ,ny | 
ve  =  ( q  Qs  i  n  <5n  y  +  vtnx}/ I  nx’nyi 
we  =  q  qCOS  5 

In  the  above,  q  and  6  are  the  velocity  magnitude  and  flow 
direction.  vt  is  selected  by  noting  the  relative  sizes 

of  9  and  the  Riemann  slip  line  direction.  The  following 

al  gori  thm  i  s  used: 


'i+l.j 


'i  +  1/2.J 


vt  i  i 


i  f  9  > 


if  9  <  Sf 
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IF? 

l  >■ 


11. 


12. 


Using  p  p  ,  u  ,  v  ,  w  ,  for  each  of  the  four  cell  edges,  fluxes 
are  evaluated  and  the  flow  field  is  advanced  using  Eqs.  (28). 

The  computational  step  is  completed  by  calculating  cell  end  areas  at 
2  +  A z  using  Eqs.  (27a)  and  decoding  to  determine  the  primitive 
variables  as  follows: 

“l.j  *  (u3/ul>U 
K/Ui)*  ,• 


Vi  ,j 


Wi  ,j 


4  1  i  ,j 


tftt 


[|i-^x|:1/2 


(36) 


where:  x  *  [(2HQu^  -  -  u2)/u^]. ^ 

Pi,j  -  <u2  - 
°1,j  *  ult 

•  9  J 

In  the  above,  p  is  positive  provided  that  x  >  1* 

The  step  size  Ac  is  chosen  using  the  CFL  condition,  which  requires  that 
the  following  condition  be  satisfied  at  every  point: 


AC 


(w2  -  a2) An 

Sc1+c2+[a2{c362+c4+c55}]i/2 


(37) 


z  — 

where:  c^  *  a  C2  -  wU 


c2  3  a  nz  -  wV 


77*2 


c,  -  (w£  -  U)"  +  (w2  -  a2)(Cv2+Cu2) 


x  y 


c4  3  (wnz-V) 2  +  (w2-a2)  (nx2+ny2' 


c5  *  2|  (w?z-U)  (wnz-V)  +  (w  -a^)  ( nxCx+nyCy, 
5  »  (An/ AC) 


NSWC  TR  86-506 


t? 


m 


*  <  »>;•> 


iV.V 


The  second  order  Godunov  scheme  is  modified  when  the  approximate  Riemann 
problem  is  used.  Steps  7  and  8  are  removed  while  the  fluxes,  determined  in 
steps  10  and  11,  are  computed  by  the  three-dimensional  version  of  Eq.  (18): 


F*  = 


a"a+(U* -If)  |  (nx,n  ) 


+  a  F  -a  F 


(a+-a‘) 


(38) 


Here,  F  ,  F",  F*  are  the  total  flux  through  the  control  volume  edges  (i.e., 
the  fluxes  of  Eq.  18  multiplied  by  the  edge  areas).  However,  U+,  IT  are  still 
the  end-area  flux  per  unit  area. 


3.3  BOUNDARY  CONDITIONS 


Five  types  of  boundary  conditions  can  be  applied  along  zone  edges, 
are  listed  in  Table  1  along  with  the  edges  on  which  each  can  occur. 


These 


The  INTERIOR  boundary  condition  allows  interior  cells  of  adjacent  zones 
to  interface.  The  second  order  Godunov  scheme  outlined  in  the  previous 
section  is  applied.  Here  it  is  presumed  that  the  mesh  is  continuous  across 
zone  boundaries.  Failure  to  use  a  continuous  mesh  will  result  in  loss  of 
second  order  accuracy. 


The  SYMMETRY  boundary  condition  permits  the  modeling  of  symmetry  planes 
and  consists  of  a  second  order  Godunov  scheme  operating  on  an  automatically 
constructed  mirror  image  state. 


Ambient  condi dtions  are  automatically  applied  to  boundary  cells  at  which 
the  FREESTREAM  condition  is  invoked.  This  procedure  advances  the  solution  as 
if  a  surface  were  adjacent  to  the  boundary  cell.  However,  at  the  end  of  the 
computational  step,  computed  properties  are  over-written  by  ambient  ones. 


The  SURFACE  boundary  conditions  allow  a  surface  to  be  simulated.  At 
cells  adjacent  to  a  surface,  the  second  order  Godunov  scheme  must  be 
modified.  The  derivatives  normal  to  the  surface  are  calculated  using  a  one¬ 
sided  difference.  Also,  the  flux  passing  through  the  cell  edge  adjacent  to 
the  surface  must  be  computed  to  reflect  the  tangent  flow  boundary  condition. 


The  slopes  normal  to  the  surface  are  determined  using  the  following  one¬ 
sided  difference  limiter: 


ap 

IT 


0  1f  < 


otherwi se 


(39) 


F'MIN 


AC 


[ - 1 . 5p i , j  +  2p2>j  -  -5P3j,  K'(p3>j  -  p2J 


K'(P2,j  -Pl,j)] 


where:  F'  =  sign(p3j  -  p2j) 
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In  smooth  regions  of  the  flow  field,  K'  is  set  to  2  while  near  fin  leading 
edge,  a  K'  value  of  zero  is  used. 

The  flux  on  the  cell  edge  adjacent  to  the  wall  is  evaluated  using  the 
procedure: 


1. 


2. 


3. 


Following  the  predictor  step,  properties  at  cell  edges  adjacent  to 
the  walls  are  computed  by  extrapolating  predicted  p,  p,  u,  v,  w  to 
the  wall  using  the  cell  slope  normal  to  the  wall. 

The  velocity  vector  determined  in  step  1  is  turned  through  either  a 
shock  or  expansion  in  plane  R.  The  appropriate  turn  aligns  the 
velocity  with  the  wall. 


The  wal  1 


F 1  /  2 ,  j  = 
The  wal 1 


is  a  steamline  and  the  flux  at  the  wall  reduces  to: 

0 

n2p 

nxP 

nyP 

flux  is  evaluated  using  the  pressure  determined  in  step  2. 


On  edges  2  and  4,  fin  surfaces  may  form  or  disappear  as  the  solution  is 
marched  down  the  length  of  the  missile.  Here  it  is  possible  for  only  part  of 
a  call  to  be  covered  by  a  surface,  as  is  illustrated  in  Figure  9.  Such  cell 
edges  are  divided  into  an  interior  and  a  surface  part.  Separate  flux  values 
are  constructed  for  each  and  summed  tc  determine  the  total  flux  passing 
through  the  cell  edge.  The  procedure  for  calculating  the  interior  and  surface 
fluxes  follows  from  the  surface  and  interior  flux  algorithms  discussed 
above.  To  implement  tnese  schemes,  it  is  necessary  to  construct  vectors 
normal  to  the  interior  and  surface  portion  of  the  cell  edges.  The  magnitude 
of  tne  surface  normals  is  equal  to  the  areas  of  the  interior  and  surface 
portions  of  the  cell  edge,  respectively.  An  estimate  of  the  areas  of  the 
surface  and  interior  portions  of  the  cell  edge  is  determined  by  dividing  the 
cell  edge  into  six  sections  as  shown  in  Figure  9.  The  total  and  surface  areas 
are  cal cul ated  from: 


tot 


■  l 
1-1 


(Ro.  -  RI>*i 


(total  area) 


(40) 


6 

l  max[(RT  -  Rr  ),  0]AZ. 


^surf 


(surface  area) 
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where:  A£. 


=  A£g  =  Az/12,  A£.  =  Az/6,  i  =  2,  3,  4,  5 


& 


.1 

,  3  ri  *TI^T(r2  -  ri> 


1 

l  A£  . 

'i,  =  r4 


Rt  =  mi n(p  (z),  R  ) 
i  o  i 

Rg  =  max(pe  (z) ,  Rj  ) 


i; 

i 


pg  (z)  =  outer  fin  tip  radial  location 
o 


pg  (z)  =  inner  fin  tip  radial  location 


$ 


r^  =  radius  at  corner  i  (see  Figure  9  for  numbering  convention). 

The  vectors  normal  to  the  surface.,  n  and  interior,  n.,  parts  of  the  cell 
edge  then  follow  from:  su  1 


surf 


n1  rtsurf  |n 
n '  I  \ot 


nI  =  n  "  nsurf 


(41a) 


(41b) 


where: 


edge  2 


[-sinx/s  +  ¥  cost,  cosx/s 


sinxi,s,  4,v]  cylindrical 


(42a) 


surf 


t-v  i  -  v 


cartesian 
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K 


i 


t 

< 


p 


edge  4 


[-sirvr/s  +  accost,  cost/s  -  sinxg  ,  oy]  cylindrical  (42b; 


surf  = 


c-v  l- 


cartesi an 


Here  v  and  o  are  evaluated  at  [zj,  +  Rg  )]  where  J  is  first  segment  for 

J  J 


which  max(RT  -  R0  ,  0)  >  e  and  R  -  R.  -  RT  +  R0  >  e,  where  e  is  a  small 
J  °J  °J  *j  “j 

number . 


The  geometry  of  edge  3  can  be  computed  to  fit  the  domain  of  dependence  of 
the  solution  when  cylindrical  coordinates  are  used.  To  invoke  this  procedure, 
the  SHOCK  boundary  conditions  are  applied.  The  geometry  of  edge  3  is 
determined  using  the  information  contained  within  the  solution  of  the  Riemann 
ProDlem  constructed  on  plane  R.  The  two  initial  states  for  this  problem  are 
tne  free-stream  and  the  abutting  cell  edge  properties.  The  cell  edge 
properties  are  obtained  by  extrapolating  the  center  properties  to  the  edge 
using  differences  normal  to  edge  3.  The  Riemann  problem  solution  features  a 


direction,  es,  which  separates  tne  free-stream  conditions  from  other  states. 


Tims  angle  marks  the  domain  of  dependence  of  the  numerical  solution  in  plane 
R.  Since  plane  R  contains  the  cell  edge  normal,  n,  the  vector  (nx,  n  )  is 
also  in  R.  A  vector  directed  along  the  9  direction  is  thus:  * 


n  = 
s 


1 


(n  *+n  2)^2 
x  y 


(nxi  +  n  j)  -  tanesk 


(43) 


Transforming  to  polar  coordinates  yields: 


h  = — — [(n  cosi>+n  sin$)e  +  (-sin^n  +cos$n  )ej  -tan9  5  (44) 

s  7777127172  x  y  '  r  '  x  y  $  s  z  v  ' 

n  +n 

'V  ' 


x  y 


Noting  that: 


*  - 

n  =  e - -  e  -  c  e 

s  r  c  $  z  z 


(45) 


and  comparing  to  the  above  yields: 


tan9  (n2  -m2  )1/2 
•  =  x  y  ' 

'z  In  cos<®  +  n  sm$j 
*  y 
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(46) 


The  procedure  for  conputing  the  geometry  of  edge  3  is  as  follows: 
Following  the  First  Step  of  Main  Algorithm 


!.  Compute  c,  using  Eq.  (46)  evaluated  with  n  and  ny  from  the  previous 
step.  Cetl  edge  properties  are  based  on  trie  celrcenter  properties  of 
the  previous  step,  and  the  differences  from  step  1  of  the  second  order 
Godunov  algorithm. 

Following  the  Predictor  (Step  4  of  Main  Algorithm) 

2.  Re-evaluate  cz  from  Eq.  (46)  using  predictor  properties. 

3.  Advance  the  snock  location  using  the  c^  value  determined  in-step  2  by 
means  of: 

cn+i  _  cn  +  (c  +c  ♦  )AZ  (47) 

2  $  4 

4.  Update  c^  values  using  a  central  difference. 

The  algoritnm  used  to  advance  the  cell  adjacent  to  tne  edge  of  the  domain 
of  dependence  is  altered  in  two  respects.  The  derivatives  in  a  direction 
normal  to  tne  edge  3  are  computed  in  a  manner  analogous  to  that  of  Eq.  39, 
witn  K'  =  2.  The  flux  passing  through  the  cell  edge  adjacent  to  the  free 
stream  is  determined  by  a  Riemann  problem  featuring  free-stream  properties  and 
aoutting  cell  edge  properties  as  the  two  initial  states. 


3.4  SEPARATION  MODELING 


At  incidences  greater  than  a  few  degrees,  the  viscous  flow  about  a 
missile  body  will  separate  and  roll  up  to  form  leeside  vortices.  Inviscid 
calculations  instead  feature  a  strong  crossflow  shock  which  may  generate 
sufficient  vorticity  to  form  leeside  vortices.  However,  the  strength  of  these 
vortices  is  much  smaller  than  that  of  measured  vortices.  Also,  the  predicted 
inviscid  pressure  distribution  differs  markedly  from  the  measured  one. 

Several  different  methods  have  been  proposed  to  model  viscous  crossflow 
separation.1*1^  The  first  of  these  alters  the  surfa.ce  velocity  direction  near 
the  estimated  separation  line,  which  is  determined  empirically.  The  second 
metnod  limits  the  crossflow  velocity  to  a  certain  maximum  value.  Both 
techniques  destroy  tne  crossflow  shock  and  produce  strong  leeside  vortices. 
However,  modeling  the  separation  line,  which  is  the  most  physically  realistic 
procedure,  yields  erratic  pressure  values  near  the  separation  line. 


The  ZEUS  code  implements  the  crossflow  limiting  approach.  The  crossflow 
plane  velocity,  u^  +  v^,  prior  to  decoding,  is  reduced  to  the  following 
maximum  Mach  number: 

,1/2/ 


M 


=  ( <*/  v  r/b) ' 


(48) 


Here  b  is  tne  local  body  radius,  a  is  the  angle  of  attack  and  r  is  the  radial 
coordinate.  The  decode  procedure,  applied  after  the  limiting  process, 
enforces  the  correct  stagnation  enthalpy  constraint.  The  above  procedure  has 
only  oeen  applied  to  turbulent  flow,  hence  the  absence  of  a  Reynolds  numDer 
deoendency . 
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SECTION  4 


RESULTS 


The  second  order  Godunov  method  has  been  applied  to  body  alone  as  well  as 
winged  configurations.  In  all  cases  the  missile  nose  was  taken  as  sharp  and 
the  initial  data  plane  was  defined  a  short  distance  from  the  nose  tip  using 
the  approximate  conical  flow  field  generator  of  Reference  9.  The  solution  was 
advanced  at  90%  of  the  CFL  step  size  and  it  was  not  necessary  to  apply  either 
special  differencing  or  artificial  viscosity.  The  constant  K.  in  the  limiter 
Eq.  31  was  set  to  unity.  Along  edges  1  and  3,  K'  of  Eq.  39  was  set  at  2, 
while  on  edges  2  and  4  it  was  set  to  0.  All  of  the  results  shown  employ  the 
approximate  Riemann  solver  of  Section  2.2.  These  results  are  nearly  identical 
to  those  obtained  using  the  approximate  Riemann  problem  of  Section  2.3. 

Figures  10  and  11  illustrate  the  calculated  surface  pressure  and 
crossflow  plane  velocities  on  a  tangent  ogive  at  a  station  6.5  calibers  from 
tne  nose  tip.  Both  clipped  and  unclipped  results  are  shown  along  with  the 
measured  surface  pressure  distribution.14  Application  of  clipping  destroys 
tne  crossflow  shock,  introduces  a  large  leeside  vortex,  and  brings  the 
calculated  leeside  surface  pressure  into  better  agreement  with  experiment.  On 
the  windward  side,  the  unclipped  surface  pressure  results  are  closer  to 
experiment. 

Calculated  surface  pressures  are  shown  in  Figure  12  for  a  cruciform  delta 
configuration  in  the  plus  roll  orientation.  Experimental  data  are  from 
Reference  15  and  were  measured  at  a  Mach  number  of  3.7,  an  incidence  of  7.8°, 
and  Reynolds  number  based  on  diameter  of  1.8(10^).  For  these  conditions, 
attached  shocks  or  expansions  occurred  at  the  fin  leading  edges.  The 
calculation  was  started  near  the  nose  tip  using  a  (18x24)  mesh  and  run  to  a 
position  slightly  forward  of  the  wing.  The  final  section  containing  the  fins 
extended  from  z  =  80  to  102  and  was  run  using  a  (36x36)  mesh.  The  calculated 
surface  pressure  agrees  well  with  experiment  over  most  of  the  fin  surface. 

The  crossflow  velocity  vectors  and  pressure  contours  at  an  axial  station  near 
the  fin  mid-cord  are  given  in  Figures  13.  Shocks  can  be  seen  attached  to  the 
fin  edges. 

Calculations  have  been  performed  on  the  two  swept  wing  configurations 
shown  in  Figure  14.  These  bodies  were  tested  at  an  incidence  of  6°  and  Mach 
numbers  of  2.5  and  4.5. ^  An  (18x18)  mesh  was  applied  forward  of  the  wing  and 
a  (36x36)  mesh  was  used  for  the  remainder  of  the  body.  Calculated  and 
measured  wing  surface  pressures  are  shown  in  Figure  14  and  agree  well  in  most 
cases.  However,  near  the  wing  leading  edge,  computed  values  are  larger  than 
measured  ones.  Figure  15  provides  measured  and  calculated  surface  pressure  on 
the  windward  side  and  leeward  side  of  the  body.  Calculated  surface  pressures 
generally  agree  well  with  experiment;  however,  discrepancies  occur  on  the  aft 
end  of  the  body  at  Mach  4.5.  On  the  windward  side,  the  pressure  rise  due  to 


the  presence  of  the  wings  is  computed  to  occur  downstream  of  the  measured  one, 
while  on  the  leeside,  predicted  pressures  exceed  experimental  values. 

However,  these  calculated  results  are  in  excellent  agreement  with  those 
computed  in  Reference  17  and  thus,  these  discrepancies  are  likely  due  to 
viscous  effects.  The  crossflow  pressure  contours  are  shown  in  Figure  16  at  an 
axial  station  upstream  of  the  wing  tip  at  Mach  2.5  and  4.5  for  the  thick  wing 
case.  A  detached  shock  is  visible  below  both  wings  and  at  the  higher  Mach 
numoer  it  is  positioned  closer  to  the  wing  surface.  This  produces  the  strong 
wing  surface  pressure  gradients  which  are  visible  at  this  Mach  number  in 
Figure  14.  At  Mach  4.5,  the  Mach  number  normal  to  the  leading  edge  is 
supersonic,  but  it  is  not  large  enough  to  produce  an  attached  shock  which 
could  turn  flow  onto  the  plane  of  the  lower  wing  surface. 

The  wing-body-tail  configuration  of  Reference  18  is  shown  in  Figure  17. 

It  features  a  highly  swept  wing  with  subsonic  leading  edge  normal  Mach 
number.  A  (36x36)  mesh  is  applied  over  the  complete  model  and  both  a 
deflected  and  undeflected  horizontal  tail  configuration  are  considered.  The 
computed  normal  force  coefficient  and  center  of  pressure  are  given  in  Figure 
17,  with  and  without  fin  deflection,  and  agree  well  with  experiment.  The 
computed  crossflow  field  at  axial  stations  near  the  wing  trailing  edge  and  the 
middle  of  the  tail  are  shown  in  Figure  18.  Figure  19  features  a  contour  plot 
of  total  pressure  loss  at  these  same  axial  stations. 


SECTION  5 


CONCLUDING  COMMENTS 


A  second  order  Godunov  method  for  tactical  missiles  in  supersonic  flow 
has  been  developed  and  applied  to  several  different  configurations.  This 
scheme  uses  a  multiple  zone  approach  to  generate  a  grid  for  finned  tactical 
missile  shapes.  It  is  cast  in  a  finite  volume  framework  and  fits  the  bow 
shock  using  the  information  contained  in  the  Riemann  problem.  Results  have 
been  applied  to  several  different  missiles,  with  and  without  fins,  and 
satisfactory  agreement  has  been  obtained  with  measurement.  It  was  not 
necessary  in  any  of  these  cases  to  use  special  procedures  or  artificial 
viscosity. 

The  merits  of  the  second  order  Godunov  scheme  can  be  illustrated  by 
comparing  it  to  the  MacCormack  finite  difference  methods  of  References  1  and  2 
(i.e.,  the  SWINT  and  MUSE  codes).  All  of  these  methods  use  a  similar  multiple 
zone  mesh  and  calculations  have  been  carried  out  on  many  of  the  shapes 
presented  in  this  report.  The  principal  advantage  of  the  Godunov  scheme 
discussed  in  this  report  is  robustness.  This  allows  complicated  cases  to  be 
coinputed  without  user  intervention  or  the  application  of  special  procedures. 
For  example,  computation  of  the  wi  ng-body-tai 1  configuration  shown  in  Figure 
17,  using  the  methods  of  References  1  or  2  requires  fine  tuning  of  the 
artificial  viscosity  level  as  well  as  careful  selection  of  the  fin 
differencing  options.  The  second  order  Godunov  method  handled  this 
configuration  without  user  intervention  or  the  application  of  special 
procedures.  In  the  case  of  a  tangent  ogive  at  incidence  (e.g. ,  Figure  10),  as 
the  mesh  is  refined,  the  MacCormack  schemes  will  fail  unless  artificial 
viscosity  is  added.  The  Godunov  scheme  sharply  resolves  the  crossflow  shock 
using  a  fine  mesh  and  needs  no  special  attention. 

A  disadvantage  of  the  Godunov  method  is  its  computational  cost  which  is 
on  the  order  of  50%  greater  than  the  MacCormack  methods.  This  problem  can  be 
reduced  by  using  the  approximate  Riemann  problem.  Computations  applying  the 
approximate  Riemann  problem  are  30%  faster  than  those  completed  with  the  full 
problem,  while  the  results  are  nearly  identical  to  those  achieved  with  the 
complete  Riemann  problem. 
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FIGURE  6.  APPROXIMATE  RIEMANN  PROBLEM 


FIGURE  10.  COMPUTEO  CROSSFLOW  FIELD  ON  A  TANGENT  OGIVE 
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1.268  R 


FIGURE  13.  COMPUTED  CROSSFLOW  VELOCITIES  AND  PRESSURE  CONTOURS  AT  STATION 
Z  =*  24.5R 


FIGURE  14.  CALCULATED  AND  MEASURED  WING  SURFACE  PRESSURES  ON  THE  SWEPT  WING  MODEL  OF  REFERENCE  16  AT  a  =  6° 


FIGURE  15.  CALCULATED  AND  MEASURED  BODY  PRESSURES  ON  THE  SWEPT  WING  MODEL  OF  REFERENCE  16  AT 
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FIGURE  18.  CROSSFLOW  VELOCITIES  AND  PRESSURES  ON  THE  MODEL  OF  FIGURE  17 
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TABLE  1.  TYPES  OF  BOUNDARY  CONDITIONS 


TYPE 

PURPOSE 

EDGES  ON  WHICH 
ALLOWED 

INTERIOR 

PROVIDE  FOR  INTERFACING 

OF  ADJACENT  ZONES 

2,  4 

SYMMETRY 

SIMULATE  A  SYMMETRY  PLANE 

2  OF  ZONE  1 

4  OF  ZONE  MAX* 

FkEESTREAM 

IMPOSE  FREESTREAM  CONDITIONS 

AT  END  OF  EACH  STEP 

ALL 

surface 

SIMULATE  A  SOLID  SURFACE 

ALL 

SHUCK 

TRACK  DOMAIN  OF  DEPENDENCE 

3 

*  MAX  -  LARGEST  ZONE  NUMBER 
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APPENDIX  A 
THE  CFL  CONDITION 


The  (CFL)  condition  requires  that  the  domain  of  dependence  of  the 
numerical  scheme  contains  the  domain  of  dependence  of  the  governing  Partial 
Differential  Equations,  in  this  case  the  Euler  equations.  The  domain  of 
dependence  of  the  Euler  equations  for  some  point  P  lies  within  the  Mach  cone 
with  vertex  at  P.  In  the  z  =  constant  plane  containing  P,  the  domain  of 
dependence  is  a  point.  Selecting  z  =  constant  planes  which  are  at  increasing 
distances  from  P,  the  size  of  the  domain  of  dependence  increases,  as  is 
illustrated  in  Figure  A-l.  The  CFL  condition  prohibits  step  sizes  in  which 
toe  Mach  cone  covers  any  area  outside  of  the  finite  difference  stencil.  The 
finite  difference  stencil  is  the  area  enclosed  by  lines  connecting  adjacent 
points  wnose  information  is  used  to  advance  properties  at  point  P. 

The  trace  of  the  Mach  cone  in  a  z  =  constant  plane  is  illustrated  in 
Figure  A-2a.  The  shaded  area  represents  the  finite  difference  stencil  of  the 
lower  order  Godunov  scheme.  This  stencil  increases  in  size  for  the  higher 
order  scheme.  However,  near  shocks  and  other  steep  gradients,  the  second 
order  scheme  degenerates  to  first  order,  and  the  CFL  step  size  calculation 
must  oe  based  on  the  more  restrictive  first  order  domain  of  dependence. 

Computation  of  the  maximum  allowable  step  size  is  more  easily 
accomplished  in  computational  coordinates  where  the  grid  is  rectangular,  as 
shown  in  Figure  A-2b.  As  can  be  seen  from  this  figure,  the_dimensi ons  of  the 
Mach  cone  cannot  exceed  1^,  1,,  lj,  1^,  in  the  n^,  h^,  n^  directions 

respecti vely .  Here  i  values  Sre: 


3  n, 


(± 


^,0,0) 


A$  An _ 

( AC^+An)1^ 


where:  A£  3  C  mesh  spacing 

An  =  n  mesh  spacing 

-sign  appl i ed  to  i  =  2,3 

while  the  n  vectors  are: 

=  (An,  A?,  nx) 


n2  =  (-An, AS,  n2  ) 


=  1,  2,  3,  4  (A-l) 


( A-2 ) 


O3  ~  (  -  An,  -  A£ , 


) 


A-l 
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n4  =  (At^-AS,  n4 


”ne  ;  component  of  each  vector  will  be  determined  below  by  considering  the 
definition  of  the  Mach  cone.  Satisfaction  of  these  four  constraints  assures 
that  the  CFL  condition  is  met. 


F 


Satisfying  the  CFL  condition  requires  determination  of  the  maximum 
dimension  of  the  Mach  cone  on  a  plane  c  =  i n  some  direction,  n.  As  an 
example,  consider  direction  n.  and  define  prane  S  with  normal  n.,  to  be 
tangent  to  the  Mach  cone.  Onthe  5  =  cn  plane,  the  point  of  tangency  between 
S  and  the  Mach  cone  is  point  T  with  coordinates  (K^AtQ+Cp>  K2A^Q+np’ 

as  shown  in  Figure  A-3.  Here  ( £  ,  ,  sp)  are  the  coordinates  of  point  P, 

(!<1,  K2)  are  constants  and  A^q  =  Ag  -  Cp.  The  Mach  cone  and  plane  S  are 

tangent  along  a  line  TP  with  direction:  (K^A^,  K^Atg.  A Cg ) . 

The  minimum  distance  between  the  projection  of  P  and  S  on  plane 
;  =  'n  represents  the  maximum  dimension,  d^ ,  of  the  Mach  cone  in 
direction  n, .  It  is  computed  as  follows: 


d.  = 


computed 

>•,  ,n,  ,  0) 
e  n 


r— 77177  •  [KlV  w  V 


Lnl  w+nl  ^ 

%  h 


[0  1  KgjAtg 

"[~n.  ~^+n,  2-j  1/2 
*5  ‘n 


(A-3) 


s. 


E 


Since  n1  is  perpendi cul ar  to  the  Mach  cone  and  the  vector  TP: 

( A-dg, ACg,  Adg)  lies  on  the  surface  of  the  Mach  cone.  Then: 

n,  *  (<-,4^,  K?A^)’  aAq)  =  0  =>  n.  K.  +  n1  K2  =  ^ 

0  OP  1  5  1  c 

and  substituting  Eq.  (A-4)  into  Eq.  (A-3)  yields: 


d,  = 


-\15Q 

~  7~—Tnj7 

( n  ^  +n  ^  ) 

5  n 


fo  satisfy  the  CFL  condition 


dj  <  ^1  * 


(A-4) 


( A-5) 


(A-6) 


Substituting  the  definition  of  t, ,  d,  from  Eqs.  (A-2)  3nd  (A-5)  respectively 
into  Eq.  (A-6)  yields: 
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To  complete  the  calculation  of  the  CFL  condition  based  on  direction  n., 
it  is  necessary  to  determine  values  of  n.  ,  n,  and  n,  which  are  consistent 

C  n  ? 

with  the  definition  of  the  Mach  cone.  In  (x,  y,  z)  space,  the  Mach  number 
normal  to  the  Mach  cone  surface  is  unity.  Thus  if  n,  is  to  be  normal  to  the 
Mach  cone,  it  must  satisfy  the  condition:  1 


(n1«  q)2  -  a2 | n1 | 2  =  0 


( A-8) 


This  condition  can  be  transformed  to  computational  space  (c,  n,  ?)  using  the 
chain  rule  applied  to  the  components  of 

nl  =  nlf5x  +  nl  nx 
x  C  n 

v  -  "lA  *  ni  \ 

y  4  n  J 

nI  =  nl/;z  +  nl  nz  +  nl 
z  ti  n  c 


luostituting  the  above  into  £q.  A-8  yields: 


a 2 [ n ,  7>  +  n.  7n]2  +•  a2n2  +  2n.  a2[n.  +n.  ] 
1  i  -  1  -  1^1 

-»  ^  c  c  ^  n 


[0.  U+n,  V+n.  w]~  =  0 
1 C  1  n  1  c 


(A- 9) 


where:  U  =  7c  •  q 
V  =  7n  •  q 

VC  =  C^i  +  cyi  +  C^ 

7n  =  n  •  +  n  •  +  n  ■ 
xi  'yJ  zk 

The  4  and  n  components  of  vector  n.  have  previously  been  specified 
as  An  and  AC,  respectively.  Substituting  these  values  into  the  above, 
a  quadratic  equation  for  the  5  component.  Solving  for  n,  /AC  produces: 

n!r 

=  |5(a2Cz-wU)  +  (a2nz-wV)|  +  a{(wcz-U)2  +  52(wnz-V)2 

+  (w2-a2)[52(Cx2  +  Cy2)  +  (nx2  +  ny2)J  +  25[(w?z-U)  (wnz-V) 

+  (w2-a2)(nx?x+nyCy)]}1/2  /(w2-a2) 


yields 


\V.\ 
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where :  5  =  | An/  ac| 

Here  the  positive  root  is  taken  and  the  absolute  value  of  the  quantity  outside 
the  radical  is  used  to  assure  that  the  largest  possible  magnitude  for  n^  has 
been  obtained.  C 

To  determine  the  maximum  possible  step  size,  the  above  calculation  must 
be  repeated  for  the  remaining  three  directions.  The  resulting  expression  for 
t.ne  CFL  condition  is  as  follows: 


_ _ (w2,  -  a2).aa _ 

•5c  i  +  c2  +  [  a z  {c^  <52  +  +  Cj5]^ 

where:  c^  =  |a2Sz-wl)j 

c2  =  | a2  n2-wV | 

c3  =  (w-;z-U)2  +  (w2-a2)(cx2+cy2) 
c4  =  (wnz-V)2  +  (w2-a2)(nx2+ny2) 
c-  =  Z\ (wC2-U)  (wnz-V)  +  (w2-a2)  (  nx  ?x+ny  Cy )  I 
5  =  An/  AS 


A-4 
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